function make_clim_fromROMS_new(roms_dir, out_dir, new_cname, newgrid, newS, new_timespan, new_timestep)
%**************************************************
% make_clim_fromROMS_new.m  
%
% 20 June 2012 SNG adapted from DAS make_clim_fromROMS.m
%
% This creates a NetCDF file to be used for ROMS forcing.  It specifically reads in 
% history files created by another ROMS run and makes a ROMS input climatology file
% with information on T, S, u, v, ubar, vbar and zeta. It extrapolates by nearest 
% neighbor to fill out the entire field (i.e. it does the preprocess_clim.m step too)
%
% INPUT: 1) roms_dir = directory where ocean_his files from old ROMS run are
%        2) out_dir = where you want your new climatology file to go
%        3) new_cname = the name of the output climatology file
%        4) newgrid = name of .nc gridfile for your new grid, generated by grid_maker.m
%        5) newS = structure S generated in run_maker.m
%        6) new_timespan = optional, can specify [start_time end_time] 2-element vector in 
%           MATLAB datenum format for a new time range to write out your climatology 
%        7) new_timestep = optional, though if you specify a new_timespan, you must also 
%           have this variable. it's the number of steps between reading in each history
%           file, so choose 24 if you want to write a new ocean_clm that has daily fields, etc.
%
%**************************************************
%
% NOTE: this is currently written to interpolate between the differing
% depths of the two runs so if your grids have very different depths this
% could be problematic, particularly in locations where the new grid is
% shallower than the old one. It may make more sense to update this to
% interpolate by sigma layers so keep this in mind!

%%
tic

clmname = [out_dir, new_cname]; %name of new climatology file

% make the directory out_dir if needed
if ~exist(out_dir,'dir');
    mkdir(out_dir);
end

%get run information
% assumes that "roms_dir" is something like:
% /Users/PM/Documents/Salish/runs/ptx_med_2005_1/OUT
% then the lines below find the string right before "OUT", which is the
% basename, the main identifier of the run
ind = strfind(roms_dir,'/');
basename = roms_dir(ind(end-1)+1:ind(end)-1);
disp(' '); disp(['basename = ',basename]);
% also figure out if we are extracting low-passed files
oname = roms_dir(ind(end)+1:end);
if strcmp(oname,'OUT_lp') || strcmp(oname,'OUT_lp_hanning') % NOTE the assumed naming convention
    basefilename = 'ocean_his_lp_';
else
    basefilename = 'ocean_his_';
end
%get the run definition
runi = roms_createRunDef('my_run',roms_dir,basefilename);

%get index of time stamps desired
if nargin == 5;     % if want same time span as original roms run, don't need new_timespan argument
    new_timestep = 1; %keep hourly as the history files probably are
    nstart = 1;
    nend = length(runi.his.ncn);
else
    nstart = dsearchn(runi.his.nctime', new_timespan(1));
    nend = dsearchn(runi.his.nctime', new_timespan(2));
end
new_timespani = nstart:new_timestep:nend;

%get the variable time from the chosen time span
vartime = NaN(1,length(new_timespani));
for ii=1:length(new_timespani)
    %get the in file
    nn = runi.his.ncn(new_timespani(ii));
    nnstring = ['000' num2str(nn)];
    nnstring = nnstring(end-3:end);
    infile = [roms_dir '/' basefilename nnstring '.nc'];
    vartime(ii) = nc_varget(infile,'ocean_time'); %define time vector
    %get a few things for later use
    if ii == 1;
        infile1 = infile;
        [G,S,~] = Z_get_basic_info(infile1);
    end
end
    
% Read in the new grid matrices & create the NetCDF output file
clim_netcdf_new(newgrid,clmname,newS.N);
lon_rho = nc_varget(newgrid,'lon_rho');
lat_rho = nc_varget(newgrid,'lat_rho');
h = nc_varget(newgrid,'h');
    
%% u,v,s,t

% set up list of variable names to create climatology files for
varname_list = {'salt';'temp';'u';'v'};
  
for ivv = 1:length(varname_list) % start of VARNAME LOOP
    varname = varname_list{ivv};
    %get variable attribute names etc.
    [V,AAlon,AAlat,~,AAlon_old,AAlat_old,AAmask_old,k] = ...
        clim_netcdf_fromROMS(varname,newgrid,G);
    [M,L] = size(AAlon);
    
    disp(['    - writing ',V.invarname,' to ',V.ncvarname]);
    
    %define rest of variables and attributes
    clim_netcdf_time(clmname,V,vartime);
    %add the time variable
    if V.do_addtimevar
        nc_varput(clmname, V.nctimename, vartime);
    end
    
    %initialize ubar and vbar when extracting u and v
    %also set which k and mask to use
    switch varname
        case 'u'
            ubar_mat = NaN(length(vartime),M,L);
        case 'v';
            vbar_mat = NaN(length(vartime),M,L);
    end
    
    %then the field
    %add attributes for the variable
    varstruct.Name = V.ncvarname;
    varstruct.Dimension = {V.nctimename, 's_rho', V.ncetaname, V.ncxiname};
    long_name = V.nclongname;
    units = V.ncunits;
    varstruct.Attribute = struct('Name', ...
        {'long_name','units'},'Value',{long_name,units});
    nc_addvar(clmname, varstruct);
    %get depth info from original grid
    indepth = S.s_rho;
    AN = length(indepth);
 
    for tt = 1:length(vartime) % start of TIME LOOP
        nn = runi.his.ncn(new_timespani(tt));
        nnstring = ['000' num2str(nn)];
        nnstring = nnstring(end-3:end);
        infile = [roms_dir '/' basefilename nnstring '.nc'];
        disp(['  tt = ',num2str(tt),' out of ', num2str(length(vartime))])
        
        if tt == 1
            %we make the new ROMS depths on the relevant grid
            hh = interp2(lon_rho, lat_rho, h, AAlon, AAlat);
            [z_rho,z_w] = Z_s2z(hh,0*hh,newS);
            hh_old = interp2(G.lon_rho, G.lat_rho, G.h, AAlon, AAlat);
            [z_rho_old, ~] = Z_s2z(hh_old,0*hh_old,S);
            DZ = abs(diff(z_w,1,1));
        end

        % reinitialize the storage arrays
        % AA is the empty matrix on the new ROMS grid in which to put
        % the field from this time step, but it has the vertical grid
        % of the original ROMS file
        AA = NaN(AN,M,L);
        % AAA is the empty matrix on the ROMS grid in which to put
        % the field from this time step, and it has the vertical grid
        % of the ROMS model.
        AAA = NaN(newS.N,M,L);
        
        % load the variable field at this time step
        clear A AAnew this_A
        % INTERPOLATE HORIZONTALLY TO ROMS GRIDS (fills whole domain,
        % if there is ANY data on that level)
        A = nc_varget(infile, varname); %read in variable from old roms file
        for zz = 1:AN % start of depth loop
            this_A = squeeze(A(zz,:,:));
            AAgood = this_A(AAmask_old);
            %fill in masked values with nearest neighbor
            AAnew = this_A;
            AAnew(~AAmask_old) = AAgood(k);
            %interpolate onto the new grid
            this_AA = interp2(AAlon_old, AAlat_old, AAnew, AAlon,AAlat);
            %save that depth level
            AA(zz,:,:) = this_AA;
        end
        
        % INTERPOLATE VERTICALLY TO new ROMS S-SURFACES
        for jj = 1:M % start of XY-LOOP
            for ii = 1:L
                this_AA = squeeze(AA(:,jj,ii));
                this_z_roms_old = squeeze(z_rho_old(:,jj,ii));
                this_z_roms = squeeze(z_rho(:,jj,ii));
                %this extraps to bottom for ncom data that is 0 sometimes
                %and include padding to the top up to zero depth
                this_AAA = interp1([this_z_roms_old; 0], [this_AA; this_AA(end)], this_z_roms,'nearest','extrap');
                AAA(:,jj,ii) = this_AAA;
            end
        end
        
        % figure out the ubar or vbar field at this time step, and write it
        % to a matrix
        switch varname
            case 'u'
                ubar_mat(tt,:,:) = squeeze(sum(AAA.*DZ,1)) ./ hh;
                if tt == length(vartime)
                    tempufile = [out_dir 'temporary_ubar_storage.mat'];
                    save(tempufile,'ubar_mat')
                end
            case 'v'
                vbar_mat(tt,:,:) = squeeze(sum(AAA.*DZ,1)) ./ hh;
                if tt == length(vartime)
                    tempvfile = [out_dir 'temporary_vbar_storage.mat'];
                    save(tempvfile,'vbar_mat')
                end
        end
        
        % write the field at one time level
        [nz,nx,ny]=size(AAA);
        nc_varput(clmname, V.ncvarname, AAA, [tt-1 0 0 0], [1 nz nx ny]);
        clear nz nx ny AAA %new_data
        
    end % end of TIME LOOP

end % end of VARNAME LOOP


%% ssh

% set up list of variable names to create climatology files for
varname = 'zeta';

%get variable attribute names etc.
[V,AAlon,AAlat,~,AAlon_old,AAlat_old,AAmask_old,k] = ...
    clim_netcdf_fromROMS(varname,newgrid,G);

disp(['    - writing ',V.invarname,' to ',V.ncvarname]);

%define rest of variables and attributes
clim_netcdf_time(clmname,V,vartime);
%add the time variable
if V.do_addtimevar
    nc_varput(clmname, V.nctimename, vartime);
end

%then the field
%add attributes for the variable
varstruct.Name = V.ncvarname;
varstruct.Dimension = {V.nctimename, V.ncetaname, V.ncxiname};
long_name = V.nclongname;
units = V.ncunits;
varstruct.Attribute = struct('Name', ...
    {'long_name','units'},'Value',{long_name,units});
nc_addvar(clmname, varstruct);

for tt = 1:length(vartime) % start of TIME LOOP
    nn = runi.his.ncn(new_timespani(tt));
    nnstring = ['000' num2str(nn)];
    nnstring = nnstring(end-3:end);
    infile = [roms_dir '/' basefilename nnstring '.nc'];
    disp(['  tt = ',num2str(tt),' out of ', num2str(length(vartime))])

    % load the variable field at this time step
    clear A AA AAnew this_A
    % INTERPOLATE HORIZONTALLY TO ROMS GRIDS (fills whole domain,
    % if there is ANY data on that level)
    A = nc_varget(infile, varname); %read in variable from old roms file
    this_A = A;
    AAgood = this_A(AAmask_old);
    %fill in masked values with nearest neighbor
    AAnew = this_A;
    AAnew(~AAmask_old) = AAgood(k);
    %interpolate onto the new grid
    this_AA = interp2(AAlon_old, AAlat_old, AAnew, AAlon,AAlat);
    %save that depth level
    AA = this_AA;

    % write the field at one time level
    [nx,ny]=size(AA);
    nc_varput(clmname, V.ncvarname, AA, [tt-1 0 0], [1 nx ny]);
    clear nx ny AA %new_data
    
end % end of TIME LOOP

%% ubar, vbar

varname_list={'ubar';'vbar'};
for ivv = 1:length(varname_list) % start of VARNAME LOOP
    varname = varname_list{ivv};    
    %get variable attribute names etc.
    [V,~,~,~,~,~,~,~] = clim_netcdf_fromROMS(varname,newgrid,G);
    
    disp(['    - writing ',V.invarname,' to ',V.ncvarname]);
    
    %define rest of variables and attributes
    clim_netcdf_time(clmname,V,vartime);
    %add the time variable
    if V.do_addtimevar
        nc_varput(clmname, V.nctimename, vartime);
    end
    
    %then the field
    %add attributes for the variable
    varstruct.Name = V.ncvarname;
    varstruct.Dimension = {V.nctimename, V.ncetaname, V.ncxiname};
    long_name = V.nclongname;
    units = V.ncunits;
    varstruct.Attribute = struct('Name', ...
        {'long_name','units'},'Value',{long_name,units});
    nc_addvar(clmname, varstruct);
    
    switch varname
        case 'ubar'
            load(tempufile)
            nc_varput(clmname,V.ncvarname,ubar_mat);
            delete(tempufile);
        case 'vbar'
            load(tempvfile)
            nc_varput(clmname,V.ncvarname,vbar_mat);
            delete(tempvfile);
    end
    
end

disp('DONE')

toc

